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Linear-response time-dependent (TD) density-functional theory (DFT) has been implemented in 
the pseudopotential wavelet-based electronic structure program BigDFT and results are compared 
against those obtained with the all-electron Gaussian-type orbital program deMon2k for the cal- 
culation of electronic absorption spectra of N2 using the TD local density approximation (LDA). 
The two programs give comparable excitation energies and absorption spectra once suitably exten- 
sive basis sets are used. Convergence of LDA density orbitals and orbital energies to the basis-set 
limit is significantly faster for BigDFT than for deMon2k. However the number of virtual orbitals 
used in TD-DFT calculations is a parameter in BigDFT, while all virtual orbitals are included in 
TD-DFT calculations in deMon2k. As a reality check, we report the x-ray crystal structure and 
the measured and calculated absorption spectrum (excitation energies and oscillator strengths) of 
the small organic molecule A'-cyclohexyl-2-(4-methoxyphenyl)imidazo[l,2-a]pyridin-3-amine. 



I. INTRODUCTION 

The last century witnessed the birth, growth, and increasing awareness of the importance of quantum mechanics 
for describing the behavior of electrons in atoms, molecules, and solids. One key to unlocking the door to widespread 
applications was the advent of scientific computing in the 1960s. New computational methods had to be developed 
— and continue to be developed — to make use of the increasing computational power. Wavclct-based methods are 
a comparative newcomer to the world of electronic structure algorithms (Refs. [l|, Q provide useful, if dated, reviews) 
but offer increasing accuracy for grid-based density-functional theory (DFT) calculations of molecular properties. 
Although DFT is a theory of the ground stationary state, excited states may be treated by the complementary time- 
dependent (TD) theory. This paper reports the first implementation of TD-DFT in a wavelet-based code — namely 
our implementation of linear response TD-DFT in the BigDFT Q program. TD-DFT results obtained with our 
implementation in BigDFT and with the Gaussian-type orbital (GTO) based program DeMon2k Q are compared 
and BigDFT is used to calculate and analyse the absorption spectrum of a recently-synthesized small organic molecule 

Modern computing requires discretization and this has been historically treated in DFT and non-DFT ab initio 
@ electronic structure calculations by basis set expansions. Because the basis set is necessarily finite, electronic 
structure algorithms for solids and molecules have developed as refinements on physically-sensible zero-order systems, 
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based upon the reasonable idea that fewer basis functions should be needed when the underlying physics is already 
approximately correct. Thus computations on solids used plane- wave basis sets partly because of analogies with the 
free-electron model of metals and partly for other reasons (ease of integral evaluation and Bloch's theorem taking 
explicit account of crystal symmetry.) Similarly computations on molecules are typically based upon the linear 
combination of atomic orbital (LCAO) approximation to molecular orbitals, with split-valence, polarization, and 
diffuse functions added as needed to go beyond this crude first approximation. A key feature of finite-basis-set 
calculations has been the use of the variational principle in order to guarantee smooth convergence from above to the 
true ground-state energy. 

Meanwhile a different set of methods developed within the larger context of engineering applications, where model- 
independent methods are important for treating a variety of complex systems. These methods involve spatial dis- 
cretization over a grid of points and include finite-difference methods, finite-element methods, and (more recently) 
wavelet methods. Use of these methods in solid-state electronic structure calculations is particularly natural in ar- 
eas such as materials science where physics meets engineering Q, in real-time applications Q and also in chemical 
engineering Q. 

Although one of the first applications of computers to scientific problems was to solve the Schrodinger equation 
by direct numerical integration using the Cooley-Cashion-Zare approach flol - [T2j and this method soon became the 
preferred way to find atomic orbitals [l3l - [l5j . quantum chemists have been slow to accept grid-based methods. There 
have been at least three reasons for this. The first is the principle that any two quantum chemistry programs should 
give the same answer to machine precision when the same calculation is performed. Historically this principle has 
been essential for debugging and assuring the consistency of the multiple sophisticated GTO-bascd programs common 
in quantum chemistry. The second reason is the fear that numerical noise would undermine the variational principle 
and make chemical accuracy (1 kcal/mol) impossible to achieve for realistic chemical problems. The first principle has 
gradually been abandoned with the wide-spread adaption in quantum chemistry of DFT and its associated grid-based 
algorithms. The second reason is simply an obsolete fear as grid-based methods are now more precise and free of 
numerical error than ever before. A third reason for the slow acceptance of grid-based methods in quantum chemistry 
is the conviction that the continuum-like unoccupied orbitals typically produced by grid-based methods are not the 
most efficient orbital-basis set to use when large configuration interaction (CI) or many-body perturbation theory 
(MBPT) type expansions are needed for the accurate description of electron correlation. Rather it is more useful to 
replace continuum-like unoccupied orbitals with natural orbitals or with other orbitals localized in the same region 
of space as the occupied orbitals since this is where electron correlation occurs. We discuss this problem further in 
Sec. El 

Wavelet-based codes represent some of the latest evolutions in grid-based methods for the problem of electronic 
structure calculations. Daubechies classic 1992 book [l6| helped to establish wavelet theory as the powerful, flexible, 
and still-growing toolbox we know today. Wavelet theory has deep roots in applied mathematics and computational 
theory, and has benefited from the various points of view and expertises of researchers in markedly different but 
complementary disciplines. (Ingrid Daubechies describes this aspect of wavelet theory rather well in a 1996 view- 
point article [131 •) Wavelet theory somewhat resembles the older Fourier theory in its use of continuous and discrete 
transforms and underlying grids. However wavelet theory is designed to embody the powerful idea of multiresolution 
analysis (MRA) that coarse features are large-scale objects while fine-scale features tend to be more localized. Compu- 
tationally, a low resolution description is afforded by a set of so-called "scaling functions" placed at nodes of a coarse 
spatial grid. Note that scaling functions are typically only placed where they are needed on this grid. Wavelets or 
"detail functions" are then added adaptively on a finer grid in regions of space where more resolution is needed. This 
provides yet another attractive feature for grid-based calculations, namely flexible boundary conditions. A subtler, 
and also highly desirable feature, is that the coefficients of the scaling and wavelet expansions are of comparable size 
and so properly reflect an even balance between different length scales. The reader interested in more information 
about wavelets is referred to the classic book of Daubechies [l6| and for applications in theoretical physics is referred 
to Refs. 

The first applications of wavelet theory to solving the Schrodinger equation began in the mid-1990s 0. (Interest 
was also expressed by quantum chemists at around the same time fl9l. l20j|.) At least two important wavelet-based 
codes have been developed for solving the Kohn-Sham equation of the traditional Hohcnbcrg-Kohn-Sham ground-state 




OMe 



FIG. 1: A f -cyclohexyl-2-(4-methoxyphenyl)imidazo[l,2-a]pyridin-3-amine (Flugi 6). 

DFT UlEi]. These are Madness [H[ and BigDFT 0. This paper concerns BigDFT. 

Since Hohenberg-Kohn-Sham DFT is a ground-state theory, a different theory is needed to treat electronic excited 



states and, in particular, to calculate absorption spectra. The TD-DFT formalism [24[ complements that of DFT by 
laying the ground work for calculating the time-dependent response of the charge density to an applied perturbation. 
Excitation spectra may then be calculated through linear response theory using, for example, the equations developed 
by one of us [25} . Implementing these equations in a wavelet code is not entirely straightforward since integral 
cvalulation is performed differently than in traditional TD-DFT codes. Thus a key point in the present paper is how 
we handle integral evaluation in our implementation of TD-DFT in BigDFT. 

Our implementation is first validated by comparison against TD-DFT calculations with the GTO-based program 
deMon2k for the historically-important test case of N2 (2g,|27|. This allows us to see and discuss some of the pros and 
cons of wavelets versus GTOs. We then go on to apply the method to a real-world application, namely the calculation 
of the absorption spectrum of a molecule of potential interest as a fluorescent probe in biological applications. This 
molecule, A r -cyclohcxyl-2-(4-methoxyphcnyl)imidazo[l,2-a]pyridin-3-amine (Fig. [I}, will simply be referred to as Flugi 
6 Q. Since this molecule has not been thoroughly characterized before we have included an experimental section 
describing our determination of its x-ray crystal structure. 

This paper is organized as follows: The basic equations of TD-DFT are reviewed in the next section. In Sec. IIII1 
we briefly review the idea of wavelets, explain how we handle the integral evaluation in our implementation of TD- 
DFT, and Section HVl discuss the pros and cons of the wavelet implementation in the context of N2 calculations with 
BigDFT and the GTO-based code DeMon2k. Section [V] reports the x-ray crystal structure, and experimental and 
calculated UV/visible absorption spectrum of Flugi 6. An assignment is also made of the peaks appearing in the 
spectrum. The final section contains our concluding discussion. 



II. TIME-DEPENDENT DENSITY-FUNCTIONAL THEORY 

As evidenced by the numerous review articles [§!■ I25I |281 - |43T | and books [43 - l46j written on the subject, TD-DFT is 
now such a well-established formalism that little additional review seems necessary. Nevertheless the purpose of this 
section is to recall a few key equations in order to keep the present paper reasonably self-contained and to introduce 
notation. 

Time-dependent DFT builds upon the Kohn-Sham formulation of ground-state DFT [22| . In its modern formulation, 
Kohn-Sham DFT is spin-DFT so that the total energy involves an exchange-correlation (xc) energy which depends 
upon the density p a of spin-a electrons and the density p$ of spin-/? electrons. These are obtained as the sum of the 
densities of the occupied orbitals of each spin, 

/V(r) = eJ2 n pM P *(r)\ 2 1 (2-1) 

p 

where n p<T is the occupation number of orbital ip pa and, the total charge density is given as the sum of the charge 
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densities of two spin components. The Kohn-Sham orbitals are obtained by solving the Kohn-Sham equation, 

— V 2 + v°[p a , Pf) }(r)\ W( r ) - e pCT </v(r) . (2.2) 

This involves a single-particle potential, 

<[Po, ^](r) = tw(r) + u H [p](r) + < c [p Q , ^](r) . (2.3) 

which is the sum of the external potential w ext (i.e., the interaction of the electrons with the nuclei and any applied 
electronic fields), the Hartree potential, 

»*)=/ ^ rfr '< ( 2 - 4 ) 
which may alternatively be obtained as the solution of Poisson's equation, 

V 2 v H [p](r) = -Anp(r) , (2.5) 
and the xc-potential which is just the functional derivative of the xc-energy, 

v X c\P*,Pi, W = — ^— 7-t — ■ (2.6) 

Together ^[p] and v£ c [p Q , p^] constitute the self-consistent field (SCF), v% xc [p a ,pp\- Modern DFT often uses a type 
of generalized Kohn-Sham formalism where the SCF may contain some orbital dependence due to (say) inclusion of 
some fraction of Hartree-Fock exchange and the external potential may even include a nonlocal part. In particular 
this latter generalization is the case when nonlocal pseudopotcntials are employed, as is the case in BigDFT. 

The external potential is time-dependent in TD-DFT, and the time-dependent Kohn-Sham potential has the form, 

ti 2 \ d 

—V 2 + v°[p a ,p p ;$ ; # ](r, t) ip p<7 (r, *) = ^<M r ) , (2.7) 

where the external potential, 

< [p a ,p p ;® ; tf ] (r,t) = v ext (r,t)+ v H [p] (r , t) + v a xc [p a , Pf} ; $ ; *o] (r, f ) . (2.8) 

The dependence of the xc part on the initial interacting (^o) an( l nonintcracting ($o) wave functions may be eliminated 
by using the first Hohenberg-Kohn theorem (2lj if the initial state is the ground stationary state. This is case when 
seeking the linear response of the ground stationary state to an applied electronic field. The TD charge density, 

p(r,t) = ej2n pa \ip pa {r,t)\ 2 , (2.9) 

pa 

then suffices to calculate the induced dipolc moment (at least for finite systems such as molecules) and hence the 
dynamic polarizability, 

e 2 f 
m e {ujf - uj*) 

whose poles give the excitation energies of the system, 

ha Jl = E I -E , (2.11) 
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and whose residues give the corresponding oscillator strengths 

2uj I m 



fi 



(2-12) 



The excitation energies and oscillator strengths are often presented in the form of a stick spectrum consisting of 
lines of height fj located at the associated w,. What is actually measured is the molar extinction coefficient e, in 
Beer's law. To a first approximation, this is related to the spectral function. 



(2-13) 



by, 



irN A e 2 



2eom e cln(fO) 



S (<»).. 



(2.14) 



in SI units [47H49| . Finite spectrometer resolution, vibrational, and solvent effects are traditionally approximated 
by replacing the Dirac delta functions in Eq. (|2 . 1 3[) with a Gaussian whose full-width at half-maximum (FWHM) 
is selected to match the experimental spectrum and we will do the same. Solvent effects may also shift excitation 
energies (50l - [56| and affect oscillator strengths [57i-l67|. While these effects are not necessarily small (vide infra), we 
will simply ignore them when comparing with experimental data as this is a reasonable first approximation. 
Practical TD-DFT calculations typically make use of the TD-DFT adiabatic approximation, 



;[p a ,Pf,}(r,t) = 



SExclp^p 



SpUv) 



(2-15) 



because even less is known about the exact xc-potential in TD-DFT than is the case for regular DFT. Here pl-(j) is 
p(r,t) regarded at fixed t as a function of r = (a;, y, z). If in addition, the xc-energy functional is of the local density 
approximation (LDA) or generalized gradient approximation (GGA) type, the TD-DFT adiabatic approximation is 
typically found to work well for low-energy excitations of dominant one-hole/one-particle character which are not too 
delocalized in space and do not involve too much charge transfer. Several extensions of classic TD-DFT have been 
found to be useful for going beyond these restrictions (see for example Ref. [42|, [43j for a review) . 

One of us used a density-matrix formalism to develop a random-phase approximation (RPA) like formalism for 
the calculation of absorption spectra from the dynamic polarizability [Eq. (|2.10p ] (2f|. This allowed the calculation 
of TD-DFT spectra to be quickly implemented in a wide variety of quantum chemistry programs since most of the 
available computational framework was already in place. The precise equation that needs to be solved is, within the 
TD-DFT adiabatic approximation, 



(2.16) 



A 


B " 




= fruj 


' 1 




B* 


A* 


(f) 


-1 


(?) 



Here, 



and the coupling matrix is given by, 



^(r)W(r) 



5 2 E. r 



r - r' 



5 Pa (v)8p T (v') 



(2.17) 



(2.18) 



6 



where we are making use of the index convention that orbitals a, 6, • • • , g, h are unoccupied, orbitals i,j, k, I, m, n are 
occupied, and orbitals o,p, ■ ■ ■ ,x,y,z are free to be either occupied or unoccupied. In the case of an LDA or GGA, 
Eq. (|2.16|) may be rearranged to give the lower-dimensional matrix equation, 

flF = h 2 uj 2 F, (2.19) 

where 

^iaajbr — &ia& jb^ar^aa ^-ia) ~t~ ' (^aa ^ia^-^iacrjbr \/ {^aa Qcr) • (2.20) 

Alternatively the Tamm-Dancoff approximation (TDA) is sometimes found to be useful [68| . 

AX = hu;X. (2.21) 

This is particularly the case when it is necessary to attenuate the effect of spin-instabilities on potential energy surfaces 
when investigating photoprocesses [69| . 



III. IMPLEMENTATION IN BIGDFT 



Our implementation of the equations of the previous section in BigDFT is described in this section and their 
validation by comparison against calculations with deMon2k is described in the next section. 

BigDFT solves the Kohn-Sham equation [Eq. (|2.2[) ] in the pseudopotential approximation. The main difference is 
that the external potential part, v ext , of the Kohn-Sham potential, vj, is manipulated so as to smooth the behavior 
of the Kohn-Sham orbitals in the core region near the nuclei while preserving the form of the Kohn-Sham orbitals 
outside the core region. This is done through the use of Goedecker-Teter-Hutter (GTH) pseudopotentials [z3| in order 
to avoid "wasting" wavelets on describing the nuclear cusp. These include both a local and nonlocal part and are used 
for all atoms (even hydrogen). Several different functionals are available for the xc-energy. However we will only be 
conside ring the LDA functional with Teter's Pade approximation [7(| of Ceperley and Alder's quantum Monte Carlo 
results [71[ since the xc-kernel, f xc , is thus far only implemented in BigDFT at the LDA level. 

There are two fundamental functions in Daubechies family: the scaling function <fi(x) and the wavelet <p(x) (see Fig. 
[2j) Note that both types of function are localized with compact support. The full basis set can be obtained from all 
translations by a certain grid spacing h of the scaling and wavelet functions centered at the origin. These functions 
satisfy the fundamental defining equations, 

II! 

cf>(x) = V2 J2 M(2z-j), 
j=l-m 
m 

<p{x) = V2 9j<l>(2x-j). (3.1) 

3=1— m 

which relate the basis functions on a grid with spacing h and another one with spacing h/2. The coefficients, hj 
and gj, consitute the so-called "filters" which define the wavelet family of order m. These coefficients satisfy the 
relations, . hj = 1 and gj = (— l) J /i_j+i. Equation (|3.ip is very important since it means that a scaling-function 
basis defined over a fine grid of spacing h/2 may be replaced by combining a scaling-function basis over a coarse grid 
of spacing h with a wavelet basis defined over the fine grid of spacing h/2. This then gives us the liberty to begin 
with a coarse description in terms of scaling functions and then add wavelets only where a more refined description 
is needed. In principle the refined wavelet description may be further refined by adding higher-order wavelets where 
needed. However in BigDFT we restricted ourselves to just two levels: coarse and fine associated respectively with 
scaling functions and wavelets. 
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FIG. 2: Daubechies scaling function (j>(x) and wavelet tp(x) of order 16. 



For a three-dimensional description, the simplest basis set is obtained by a tensor product of one-dimensional 
basis functions. For a two resolution level description, the coarse degrees of freedom are expanded by a single three 
dimensional function, 0^ 4 4 (r) , while the fine degrees of freedom can be expressed by adding another seven basis 
functions, cjf^ - 2 j 3 (r), which include tensor products with one-dimensional wavelet functions. Thus, the Kohn-Sham 
wave function ^(r) is of the form 

*W= E <&.«,.«,<.«..«,(*) + E E4^3^3«- (3-2) 

ii,i2,i3 31,32,33 v\ 

The sum over £1,^2,23 runs over all the grid points contained in the low- resolution regions and the sum over J2J3 
runs over all the points contained in the (generally smaller) high resolution regions. Each wave function is then 
described by a set of coefficients {c^ - 2 J3 }, v = 0, ...,7. Only the nonzero scaling function and wavelet coefficients 
are stored. The data is thus compressed. The basis set being orthogonal, several operations such as scalar products 
among different orbitals and between orbitals and the projectors of the nonlocal pseudopotential can be directly 
carried out in this compressed form. In addition to raw Daubechies scaling functions, practical applications make use 
of autocorrelated functions to make interpolating scaling functions (ISF) [z2]. In particular, as shown by A. Neelov 
and S. Goedecker (73|, the local potential matrix elements approximated using the linear combination of such ISFs 
is exact for polynomial expansions up to 7th order and the corresponding Kohn-Sham densities can be calculated by 
the real-valued coefficients on the grid points. 

Although not every grid point is associated with a basis function and the fine grid is only used in some regions 
of space, the Daubechies basis set is still very large. This means that full diagonalization of the Kohn-Sham orbital 
Hamiltonian is not possible. Instead the direct minimisation method [zl, [75| is used to obtain the occupied orbitals. 
This is in contrast with the GTO-bascd deMon2k program which will be described below in which full diagonalization 
of the Kohn-Sham orbital Hamiltonian matrix is carried out within the finite GTO basis set. 

A key point to review because of its importance for our implementation of TD-DFT in BigDFT is the Poisson 
solver used to treat the Hartree part of the potential, v H . Although this Poisson solver has been discussed elsewhere 
[Zfil IZZl ; we briefly review how it works in order to keep this article somewhat self-contained. The Hartree potential 
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is evaluated as, 

v H (r) = J G(\r-r'\)p(v')dv', (3.3) 

where G(r) = l/r is the Green's function for the Poisson equation, namely just the Coulomb potential. The density 
and potential are expanded in a set of interpolating scaling functions, 

P( r ) = ^i( x )^2(y)<M 2 )Pn,i2,*3 

h, 12,13 

" H (r) = 4) il (x)4>. l2 (y)4> l3 {z)v lul2 ^ 3 , (3.4) 

associated with the same grid of points, rj li j 2i j 3 , in real space. In particular, the charge density coefficients, Pi lt i 2 ,i 3 = 
p( r ii,i 2 ,i 3 )- Then, 

v iiM,i3 = y ] Gj 1 —i lt j 2 —i 3> j 3 —i 3 pj lt j 2t j 3 , (3-5) 

where the quantity, 

J l r *l,*2,*3 1 I 

is translationally invariant by construction. Since Eq. (|3.5p has the form of a three dimensional convolution, it may be 
efficiently evaluated by using appropriate parallelized fast Fourier transform algorithms at the cost of only 0(N In N) 
operations. The calculation of matrix elements of the Green's function G(r) = l/r is simplified by using a separable 
approximation in terms of Gaussians, 



i«£e-^ 2 Cfc , (3.7) 



r 

k 



so that all the complicated 3-dimensional integrals are reduced to products of simpler 1-dimensional integrals. For 
more information about BigDFT, the reader is referred to the program website [1[ and to various publications [7614811 ]. 

We are now in a position to understand the construction of the coupling matrix in our implementation of TD-DFT 
in BigDFT, which we split into the Hartree and exchange-correlation parts, 

Kaia.bjr = ^aia.bjr ^aja,bjr ■ (3-8) 

Instead of calculating the Hartree part of coupling matrix directly as, 

4W=/ |v^(r)^ CT (r)^^Vv(r')^ T (rVrdr', (3.9) 
we express the coupling matrix element as, 

KioM-r = J C('^W%W*> (3-10) 

where, 

<w(r)= / f^dr', (3.11) 
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and, 

p aiCT (r)=C CT (rhMr). (3.12) 

The advantage of doing this is that, although p a io and v a ia are neither real physical charge densities nor real physical 
potentials, they still satisfy the Poisson equation, 

V 2 ?w(r) = -47rp aiCT (r), (3.13) 

and we can make use of whichever of the efficient wavelet-based Poisson solvers already available in BigDFT, is 
appropriate for the boundary conditions of our physical problem. 

Once the solution of Poisson's equation, w aj(T (r), is known, we can then calculate the Hartree part of the kernel 
according to Eq. (|3.10[) . Inclusion of the exchange-correlation kernel is accomplished by evaluating, 

K a ia,bjT = J M ai v(r)p bjT (r)dr, (3.14) 

where, 

M aia {r) = v ma (v) + [ p ailT {r')f^ T {r, r') dr' . (3.15) 



We note that /^ T (r, r') = f^'J{v, r')S(r — r') for the LDA, so that no integral need actually be carried out in evaluating 
Maicrir). The integral in Eq. (|3.14p is, of course, carried out numerically in practice as a discrete. 



IV. VALIDATION 



Having implemented TD-DFT in BigDFT we now desire to validate our implementation by testing it against 
another program in which TD-DFT is already implemented, namely the all-electron GTO-based program deMon2k 
DeMon2k resembles a typical GTO-based quantum chemistry program in that all the integrals other than the 
xc-integrals, can be evaluated analytically. In particular, deMon2k has the important advantage that it accepts the 
popular GTO basis sets common in quantum chemistry and so can benefit from the experience in basis set construction 
of a large community built up over the past 50 years or so. In the following, we have chosen to use the well-known 
correlation-consistent basis sets for this study [HI, [&|. (Note, however, that the correlation-consistent basis sets used 
in DeMon2k lack / and g functions but are otherwise exactly the same as the usual ones.) The advantage of using 
these particular basis sets is that there is a clear hierarchy as to quality. 

An exception to the rule that integrals are evaluated analytically in deMon2k are the xc-integrals (for the xc-energy, 
xc-potential, and xc-kernel) which are evaluated numerically over a Becke atom-centered grid. This is important be- 
cause the relative simplicity of evaluating integrals over a grid has allowed the rapid implemenation of new functionals 
as they were introduced. We made use of the fine fixed grid in our calculations. 

As described so far, deMon2k should have 0(N 4 ) scaling because of the need to evaluate 4-center integrals. 
Instead deMon2k uses a second atom-centered auxiliary GTO basis to expand the charge density. This allows the 
the elimination of all 4-center integrals so that only 3-center integrals remain for a formal 0(N 3 ) scaling. In practice, 
integral prescreening leads to 0(N M ) scaling where M is typically between 2 and 3. We made use of the A3 auxiliary 
basis set from the deMon2k automated auxiliary basis set library. 

All calculations were performed using standard deMon2k default criteria. Although full TD-LDA calculations are 
possible with deMon2k, the TD-LDA calculations reported here all made use of the TDA. The chosen test molecule 
was N2 with an optimized bond length of 1.115 A. This molecule was chosen partly because of its small size but also 
because of the large number of excited states which are well characterized (see Refs. [26|, [27| and references therein.) 

Unlike TD Hartree-Fock (or configuration interaction singles) calculations, TD-LDA calculations are prepreparcd 
to describe excitation processes in the sense that the occupied and unoccupied orbitals both see the same number of 
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TABLE I: Basis set dependence of the HOMO and LUMO energies and of the HOMO-LUMO gap (eV) calculated using 
deMon2k. 
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AUG-CC-PVTZ 


10.3953 


2 
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AUG-CC-PV5Z 
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2 
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7 
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AUG-CC-PCVDZ 
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2 
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7 


,9847 


AUG-CC-PCVTZ 
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2 
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7 


9957 


AUG-CC-PCVQZ 


10.3990 


2 


.4124 


7 


,9866 


AUG-CC-PCV5Z 


10.3985 


2 


.4136 


7 


,9849 



electrons (because they come from the same local potential). This means that there is often little relaxation — at 
least in small molecules — and a two orbital model (2f| SH 0, [HI often provides a good first approximation to the 
singlet (hcof^ a ) and triplet (fkof_^ a ) excitation energies, 

fv^Ua = ta- et + (ia\2f H + f%° + f«f\ai) . (4.1) 
Consideration of typical integral signs and magnitudes then implies that 

< £a - < tUjlUa , ( 4 - 2 ) 

with the singlet-triplet splitting going to zero for Rydberg states (in which case the electron repulsion integrals become 
negligible due to the diffuse nature of the target orbital ip a .) 

Since orbital energy differences provide an important first estimate of TD-DFT excitation energies, we wished to 
see how rapidly they converged for BigDFT and deMon2k as the quality of the basis set was improved. Tables [I] 
and [n] show how the highest occupied molecular orbital (HOMO) - lowest unoccupied molecular orbital (LUMO) 
energy gap (Ae 

HOMO-LUMO 

) varies for each program. 

Consider first how deMon2k calculations of Ae HOM o-LUMo evolve as the basis set is improved (Table HJ). Jamorski, 
Casida, and Salahub had earlier shown that LUMO is bound for reasonable basis sets (26[ . (Small differences between 
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TABLE II: Basis set dependence of the HOMO and LUMO energies and of the HOMO-LUMO gap (eV) calculated using 
BigDFT. 

h g 7mVn c 

Low resolution 



U.4/ 0/ o 


10.3910 


2.3815 


o.uuyo 


U.4/ I/O 


10.3964 


2.3922 


8.UU42 


0.4/8/8 


10.3971 


2.3945 


8.0027 


0.4/9/8 


10.3972 


2.3951 


8.0022 


0.4/10/8 


10.3973 


2.3953 


8.0020 


0.4/11/8 


10.3972 


2.3953 


8.0019 




High resolution 




0.3/7/8 


10.3977 


2.3932 


8.0043 


0.3/8/8 


10.3984 


2.3957 


8.0027 


0.3/9/8 


10.3985 


2.3963 


8.0022 


0.3/10/8 


10.3985 


2.3965 


8.0021 


0.3/11/8 


10.3985 


2.3965 


8.0020 



"Grid spacing of the cartesian grid in atomic units. 
'Coarse grid multiplier (crmult) 
Tine grid multiplier (frmult) 

the present calculations and those in Ref. (2(| are due to gradual improvements in the grid, auxiliary basis sets, and 
convergence criteria used in the DeMon programs). Convergence to the true HOMO-LUMO LDA gap is expected 
with systematic improvement within the series: (i) double zeta plus valence polarization (DZVP) — > triple zeta plus 
valence polarization (TZVP), (ii) augmented correlation-consistent double zeta plus polarization plus diffuse on all 
atoms (AUG-CC-PCVDZ) -> AUG-CC-PCVTZ (triple zeta) -> AUG-CC-PCVQZ (quadruple zeta) -> AUG-CC- 
PCV5Z (quintuple zeta), (hi) augmented correlation-consistent valence double zeta plus polarization plus diffuse 
(AUG-CC-PVDZ) -> AUG-CC-PVTZ -> AUG-CC-PVQZ -> AUG-CC-PV5Z, (iv) correlation-consistent double zeta 
plus polarization plus tight core (CC-PCVDZ) -> CC-PCVTZ -> CC-PCVQZ -> CC-PCV5Z, and (v) correlation- 
consistent valence double zeta plus polarization on all atoms (CC-PVDZ) -> CC-PVTZ CC-PVQZ -> CC-PV5Z. 
There is a clear tendency in the correlation-consistent basis sets to tend towards values of 10.40 eV for the HOMO 
energy, 2.42 eV for the LUMO energy, and 8.01 eV for A£ HOMO _ LUMO , with adequate convergancc already achieved 
with the CC-PVTZ basis set. 

Now let us turn to BigDFT (Table HI)) . Calculations were done for two different grid values, denoted by h g = 0.3, 
0.4. (These values are the nodes of the grid in atomic units which serve as centers for the scaling function/ wavelet 
basis.) The simulation "box" has the shape of the molecule and its size is expressed in the units of the coarse grid 
multiplier (crmult) and the hnc grid multiplier (frmult) which determines the radius for the low/high resolution sphere 
around the atom. Results using 16 unoccopicd orbitals with the different wavelet basis sets are essentially identical, 
with no significant variation in the HOMO energy, the LUMO energy, and the Ae HOMO „ LUMO value of 8.00 eV between 
the high resolution combination of 0.3/11/8 and the low resolution combination of 0.4/6/8. 

The remaining differences for the HOMO energy, LUMO energy, and Aehomo-lumo calculated by the two programs, 
deMon2k and BigDFT, is more difficult to trace. For example, it might be due to the auxiliary basis approximation 
in DeMon2k or to the use of pseudopotentials in BigDFT or perhaps to still other program differences. The important 
point is that differences are remarkably small. 

We now come to the calculation of the actual excited states of N2 and the third reason alluded to in the introduction 
that quantum chemists have been slow to adapt grid-based methods. This is the concern that the very large size of basis 
sets in grid-based methods would lead to impractically-large configuration interaction expansions. Put differently, this 
concerns the basic problem of how to handle the continuum. A correct inclusion of the continuum in the formalism of 
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FIG. 3: Singlet and triplet excitation energies for N2 calculated using BigDFT 



TABLE III: Comparison of the nine lowest excitation energies of N2 (in eV) calculated using different programs and with 
experiment. 

State BigDFT" deMon2k'' ABINFT Experiment d 



1 3 E+ 


10.71 


10.39 


9.16 


12.0 


i 3 n u 


10.58 


10.38 


10.79 


11.19 




10.29 


9.99 


10.46 


10.27 




9.37 


9.36 


9.92 


9.92 


1 3 £" 


9.36 


9.36 


9.91 


9.67 


1% 


9.35 


9.10 


9.47 


9.31 


1 3 A„ 


8.93 


8.60 


9.08 


8.88 


i 3 n s 


7.69 


7.83 


7.85 


7.75 


i 3 £+ 


8.52 


8.43 


8.16 


8.04 



"Present work (TD-LDA/TDA). 
''Present work (TD-LDA/TDA). 
c From HI (TD-LDA). 
d Taken from Ref.([§3). 
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Sec. |TI] would seem to require at least approximate integrals over a quasicontinuum of unoccupied orbitals. Clearly this 
is impractical and our method only uses the first several unoccupied orbitals. This is perhaps reasonable for the lower 
excited states, given the anticipated dominance of the two-orbital model [Eq. (|4.1I) ], but may fail to be quantitative 
when relaxation or state mixing starts to become important and the two-orbital model breaks down. This problem 
is avoided in quantum chemistry where the virtual orbitals in the TD-DFT calculation have an altogether different 
meaning: they are simply there to describe the dynamic polarizability of the ground-state charge density and need 
not describe the continuum well. This is why fewer unoccupied orbitals are normally needed in quantum chemical 
applications of TD-DFT than is the case when using, say, plane waves. In Fig. [3] the lowest few excited-states are 
calculated using BigDFT. This figure shows that the excited states calculated with different grids are essentially the 
same up to -£homo- 

After 

-f H oMOi increasing the quality of the basis set by varying the simulation box sizes leads to an 
increasing collapse of the higher excited-states. This is a reflection of the fact that the TD-LDA ionization continuum 
starts at -6 HOMO which is about 5 eV low because of the incorrect asymptotic behavior of the xc-potential [27j . Of 
course, as indicated by the estimate (|4.1[) , the excitation energy is not exactly a simple orbital energy difference, but 
they are closely related. 

With these caveats, let us see how well our implementation of TD-DFT does in BigDFT. Table IIII1 lists the 
lowest nine excited states of N2 calculated with DeMon2k and BigDFT. For comparison, Table ITU also contains full 
TD-LDA (i.e, non-TDA) excitation energies obtained from ABINIT using the Perdew-Wang 92 parameterisation 
of the LDA functional [H| along with the experimental values from the literature (86[. The slight differences which 
occur between the 1 X S^ and 1 3 S^ excitation energies in the BigDFT and ABINIT calculations are an indication 
of residual numerical errors since these two states are rigorously degenerate by symmetry when using the TD-LDA 
and TD-LDA/TDA approximations: Aside from this tiny difference, it is certainly reassuring that excitation energies 
calculated with BigDFT, deMon2k, and ABINIT are quite similar. Nevertheless differences as large as 0.3 eV or 
more are found for some states. Such differences are large enough to be potentially problematic for determining the 
ordering of near-lying states. 

Finally since a large number of excited states have been calculated, it is interesting to test the assertion that 
the oscillator strength distribution should be approximately correct even above the TD-LDA ionization threshold at 
-£homo [13 ■ This is especially true since the transitions given in Table UTTI are all dark and we would like to see how 
the oscillator strengths of bright states compare. The high-energy spectra are compared in Fig. 0] using Eq. (|2.14l) . 
Clearly the spectra are in reasonable qualitative agreement. 

The above results show that this first implementation of BigDFT is quantitative — especially when results are 
dominated by bound-bound transitions — and our discussion may eventually suggest ways to go further towards 
improving the method. In the meantime, we have a method that can be used for moderate-size molecules and this is 
illustrated in the next section. 



V. APPLICATION 

A large series of fluorescent molecules of potential interest as biological markers [HI has recently been synthesized by 
combinatorial chemistry [f| . This is a method whereby large sets of similar reactions are conducted in parallel in arrays 
of spots on a single plate, thus dramatically increasing throughput when searching for molecules with a particular 
property — in this case fluorescence. We have chosen to calculate the absorption spectrum of one of these fluorescent 
molecules in preparation for future more in-depth theoretical studies of their fluorescence properties using BigDFT. 
This molecule, which we will simply refer to as Flugi 6 (because it is molecule 6 [f| among the fluorescent molecules 
prepared by the UGI reaction [90j) rather than by its full name of A^-cyclohcxyl-2-(4-methoxyphcnyl)imidazo[l,2- 
a]pyridin-3-aminc is shown in Fig. [TJ The synthesis and partial characterization of Flugi 6 has been described in Ref. 
|5[. However we go further here and report the experimental determination of its crystal structure. We then go on 
to compare the spectrum calculated with BigDFT with the measured spectrum and discuss the problem of peak 
assignment. 
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A. X-ray Crystal Structure 



Crystals of Flugi 6 (C20H23N3O, M — 321.42 g/mol) were obtained out of recrystallisation in ethyl acetate (EtOAc) 
as colorless needles suitable for X-ray diffraction. A 0.38 mm x 0.28 mm x 0.01 mm crystal was mounted on a glass 
fiber using grease and centered on a Bruker Enraf Nonius kappa charge-coupled device (CCD) detector working at 200 
K and at the monochromated (graphite) Mo K a radiation A = 0.71073 A. The crystal was found to be orthorhombic, 
Pna2i, a = 27.912(4) A , b = 5.876(2)1, c = 10.297(2) A, V = 1688.7(6) A 3 , Z = 4, D x = 1.264 g.cm" 3 , m = 0.080 
mm -1 . A total of 17700 reflections were collected using <fi and uj scans; 2853 independent reflections (Ri n t — 0.1557). 
The data were corrected for the Lorentz and polarization effects. The structure was solved by direct methods with 
SIR92 [9l| and refined against F by leastsquare method implemented by TeXsan (92|. C, N, and O atoms were refined 
anisotropically by the full matrix least-squares method. H atoms were set geometrically and recalculated before the 
last refinement cycle. The final R values for 1964 reflections with I > 2a (I) and 217 parameters are Rl = 0.0617, 
wR2 = 0.0657, goodness of fit (GOF) = 1.78 and for all 2854 unique reflections Rl = 0.0923 , wR2 = 0.0829, GOF 
= 1.85. The resultant crystal geometry is given in IIVI 

The data have been deposited with the Cambridge Crystallographic Data Centre (Reference No. CCDC 8200007). 
This material is available free of charge via www.ccdc.can.ac.uk/conts/retricving.html (or from the Cambridge 
Crystallographic Data Centre, 12 Union Road, Cambridge CB2 1EZ, UK. Fax: +44-1223-336033. E-mail: de- 
posit@ccdc.cam.ac.uk). 
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FIG. 5: Experimental geometry: carbon, orange; nitrogen, blue; oxygen, red; hydrogen, white. This geometry consists of two 
nearly planar entities, namely a nearly planar cyclohexane (CeHu-) ring and the rest of the molecule which rests in a plane 
perpendicular to the plane of the cyclohexane 

B. Spectrum 

The experimental UV/Vis spectrum was determined in dimethyl sulphur dioxide (DMSO) as described in the 
supplementary material of Ref. Q The first step to calculating the spectrum of Flugi 6 is to optimize the geometry. 
This was done with BigDFT using the 0.4/8/8 grid, beginning with the experimental crystal geometry shown in 
Fig. [5] Our optimized geometry is given in Table IVl The largest change 0.15 A from the initial guess is for the 
18th atom of hydrogen. It was verified that the optimized geometry is indeed a minimum by explicit calculation of 
vibrational frequencies. However the experimental geometry is not identical to our calculated gas phase geometry as 
confirmed by the presence of three imaginary vibrational frequencies calculated for the (unoptimized) experimental 
geometry. 

The TD-LDA absorption spectrum of Flugi 6 was then calculated at the optimized geometry using our new imple- 
mentation of TD-DFT in BigDFT. In addition to the previously mentioned computational details, the calculation 
used 60 unoccupied orbitals within the TDA. The excited-states were obtained using full-matrix diagonalization of 
the TD-DFT part. The theoretical spectrum was calculated using Eq. (12. 141) using a FWHM of 0.25 eV and then 
transformed to a wavelength scale using our spectrum convolution program |93f . Because we were restricted to a more 
limited number of unoccupied orbitals than in the N2 test case, there is some concern that our calculated spectrum 
might change if a larger number of unoccupied orbitals is included. However the comparison of the theoretical and 
experimental results shown in Fig. [S] is reasonable. This is especially true when it is kept in mind that we are com- 
paring gas-phase theory with an experimental spectrum measured in a polar solvent DMSO. Notice the presence of a 
larger peak at 4.6-5.0 eV, a smaller peak at 3.2-3.5 cV, and a shoulder inbetween near 4.3 eV. 

In contrast to experience with ruthenium complexes (to name but one example), Eq. (|2.14l) with gaussian broadening 
does not suffice in the present case to give good agreement between theoretical and experimental molar extinction 
cocfficincts. This is why the theoretical curve has been divided magnitude by a factor of ten. We believe that this 
discrepancy is in part due to the aforementioned solvent effects on excitation energies and oscillator strengths which 
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FIG. 6: Comparison of theoretical and measured absorption spectra for Flugi 6 (left y-axis). The magnitude of the BigDFT 
curve has been divided by a factor of ten (see text). Both theoretical and experimental curves show qualitative agreement with 
the oscillator strength stick-spectra which however is in different units (right y-axis). 

we have chosen to neglect and in part due to the possible presence of multiple conformcrs in the room temperature 
experiment. This latter hypothesis might be tested by expensive dynamics calculations, but this far beyond the scope 
of the present work. Nevertheless, even without dynamics we find this level of agreement to be encouraging and now 
go on to further analyze our calculated spectrum. The calculated stick spectrum is also shown in Fig. [6] It is now 
clear that the small peak at 3.2-3.5 cV corresponds to two transitions, that the shoulder near 4.3 cV corresponds 
to three transitions, and that the large peak at 4.6-5.0 eV corresponds to several electric excited states. Table IVII 
provides a more detailed analysis. All of the transitions are below the onset of TD-LDA ionization continuum at 
-ehomo = 4.8713 cV, which is artificially low compared to the true ionization energy [27| . As mentioned in Sec. IIII1 
unlike TD Hartree-Fock (or configuration interaction singles) calculations, TD-LDA calculations arc prcprcpared to 
describe excitation processes in the sense that the occupied and unoccupied orbitals both see the same number of 
electrons (because they come from the same local potential.) This means that there is often little relaxation - 
at least in small molecules — and a two orbital model (25| provides a good first approximation to the excitation 
energy [Eqs. (|4.1[) and (|4.2[) ]. The TDA configuration interaction coefficient is then determined by spin coupling and 
is given by l/v2 = 0.707. Tabic fVj shows significant deviations from this theoretical value, suggesting significant 
relaxation effects may be taking place. Visualization of the HOMO and LUMO suggests that relaxation is important 
here and might help to explain why the first peak is at slightly too low an energy compared to the first experimental 
peak in the absorption spectrum |8J]. Nevertheless, the energy of HOMO — > LUMO dominated singlet transition at 
3.21 eV exceeds the simple difference of HOMO and LUMO molecular orbital energies of 2.80 eV as expected from 
the domination of the Hartree term (ia|/#|ai) over the two xc terms. Two interesting features, which will not be 
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pursued in the present paper, are the absence of oscillator strength for the 1 {H 1 L + 1) transition and the indication 
of significant configuration mixing for the first and seventh transitions which both borrow from 1 (H,L) and for the 
third and fifth transitions which both borrow from {H, L + 2). 



VI. CONCLUSION 



Grid-based methods have long been regarded with skepticism by the quantum chemistry community, but have 
now been accepted in the form of the grids used to evaluate xc-integrals in the DFT part of most quantum chemistry 
codes. We believe that an even greater acceptance of grid-based methods may prove useful as theoretical solid-state and 
chemical physics strive to meet on "neutral ground" at the nanometer scale. Acceptance will be aided by continuing 
advances in grid-based methods with wavelets being of particular interest here. At the same time, new features should 
be added to grid-based electronic structure codes in order to make them more useful for chemical applications. This 
paper represents a small step in that direction. In particular, we have presented the first implementation of TD-DFT 
in a wavelet-based code, namely in BigDFT. 

While BigDFT is designed for routine calculations on systems containing many hundreds of atoms, this first 
implementation of TD-DFT in BigDFT is not yet ready for these more ambitious applications. Rather, we wished to 
bring out the pros and cons of wavelet-based TD-DFT by comparing against a GTO-based quantum chemistry code 
(in this case, deMon2k) and by an example application showing how our implementation in BigDFT can be useful 
in analyzing the spectrum of a molecule of contemporary experimental interest. 

A factor in favor of the wavelet-based approach is the rapidity of convergence of the bound orbitals and orbital 
energies with respect to refinements of the wavelet basis and associated grid. While orbital results from the all-electron 
deMon2k code are not (and should not) be in exact agreement with those of the pseudopotential BigDFT code, 
the results are really quite close when sufficiently large basis sets are used. In the case of the GTO-based code, tight 
basis set convergence typically requires going beyond at least the triple-zeta-valence-plus-polarization (TZVP) level. 
In contrast, adequate orbital convergence is easily obtained with BigDFT using the default wavelet basis set and 
grid, with further refinements leading to only minor improvement. This comes close to the quantum chemists' dream 
of calculations free of errors due to basis set incompleteness. 

Interestingly, problems which could be envisioned with this first implementation of TD-DFT in BigDFT either 
did not arise or did not seem to be serious. The worry was that the unbound orbitals of the molecule arc continuum 
orbitals whose description is apparently only limited by the boundaries of the box defined by the coarse grid. In 
principle, for an infinitely large box, there are an infinite number of unoccupied orbitals in even a small energy band 
and all of these would seem to need to be taken into account even for describing transitions below the TD-DFT 
ionization threshold at -£homo- This is a doubly large worry because the number of unoccupied orbitals to be 
calculated is limited by an input parameter, making a double convergence problem (number of virtuals and box size.) 
Nevertheless our calculations show that the implementation in BigDFT works correctly, giving quite reasonable 
results when compared with deMon2k and with experimental results for transitions below -chomo- There are at 
least two probable reasons that the anticipated problems are not seen here. The first is the tendency, at least for small 
molecules, for excitations to be dominated by bound-to-bound transitions involving two or only a few orbitals. This is 
especially true for the LDA and GGA, but will gradually breakdown with the inclusion of Hartree-Fock exchange where 
orbital relaxation becomes more important. The second reason that we may not see the expected problems associated 
with continuum-type unoccupied orbitals is that putting the molecule in a box acts much like atom-centered GTOs 
in the sense that it keeps wavelet basis functions close to those regions of space where electron density is high and so 
can be efficiently used to describe the dynamic response of the charge density, whose description is the fundamental 
key to extracting spectra in TD-DFT. Nevertheless it should be mentioned that there are alternative algorithms in 
TD-DFT such as the modified Sternheimer equation and the Green's function approach which avoid explicit reference 
to unoccupied orbitals 44[. These may be worth exploring in future development work of wavelet-based TD-DFT. 
However our first priority will be to implement analytic derivatives for TD-DFT excited states which arc needed in 
modeling fluorescence. 

It should also be mentioned that implementing TD-DFT is a step along the way to implementing MBPT methods 
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from solid-state physics, namely the GW one-particle Green's function approximation and the Bethe-Salpeter equation 
approach to the two-particle Green's functio n 19411 ■ Such work is already in progress in BigDFT. Of course, the 
expected collapse of the TD-DFT continuum |27| was seen above -chomo (Fig. EJ with increasing box size, though 
the spectrum remains qualitatively correct (Fig. @|. 

The application to Flugi 6 presented in this paper provides a concrete reality check on the usefulness of TD-DFT for 
practical applications. Geometry optimization was simplified by beginning with an x-ray structure, but solvent and 
dynamics effects were ignored. All in all the final result may be qualified as semiquantitative but useful. In fact, we 
have also carried out preliminary calculations of absorption spectra for five other members of the Flugi combinatorial 
chemistry series (not reported here.) For these molecules, trends in the energies of the first experimental absorption 
peaks do parallel trends in the calculated first absorption peak as well as the HOMO-LUMO energy difference. 
However, in the absence of x-ray crystal data, the large number of possible molecular configurations merits further 
exploration, especially since the LDA may not correctly order these states. Hydrogen bonding with the solvent should 
also be considered in some cases. For these reasons it seems best to reserve the calculation of these spectra, comparison 
with experimental spectra, and assignment of transitions for a future paper. 
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TABLE IV: Experimental geomentry (Cartesian coordinates in A) for the Flugi 6 
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TABLE V: DFT optimized geometries (Cartesian coordinates in A) of Flugi 6. Calculations performed at the LDA level of 
theory. 
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TABLE VI: Singlet excitation energies (fkJi, in eV) up to — ehomo = 4.8713 eV, oscillator strength (//, unitless) and assignment. 
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a H and L stand respectively for HOMO and LUMO. 
''Configuration interaction (CI) expansion coefficient. 



